function [theta_gm, phi_gm, varargout] = gg2gm_new(theta_gg, phi_gg, varargin)
% transformation between geographic and geomagnetic (dipole) co-ordinates
%          and components
%
% theta_gm                = gg2gm(theta_gg, phi_gg)
% [theta_gm, phi_gm]      = gg2gm(theta_gg, phi_gg)
% [theta_out, phi_out]      = gg2gm(theta_in, phi_in, i_trans)
%                           i_trans = +1: gg -> gm
%                           i_trans = -1: gm -> gg
% [theta, phi, B_theta, B_phi]  = gg2gm(theta_gm, phi_gm, i_trans, B_theta, B_phi)
%
% theta(:) co-latitude [deg]
% phi(:)   longitude [deg]
%
% March 2003, Nils Olsen, DSRI

if nargout == 3; error('This is the new version of gg2gm. You are probably looking for the old one [theta_gm, phi_gm, psi] = gg2gm(theta_gg, phi_gg)'); end;

rad = pi/180;
if nargin < 3
    i_trans = +1;
else
    i_trans = varargin{1};
end
if abs(i_trans) ~=1; error('i_trans should be +1 or -1'); end

if nargin > 3
    if nargin == 5 && nargout == 4
    else
        error('5 input and 4 output variables required')
    end
end

% Initialization: coordinates of geomagnetic North pole
theta_b = 11.32;
phi_b     = 289.59;

s_p_b = sin(phi_b*rad);
c_p_b = cos(phi_b*rad);
c_t_b = cos(theta_b*rad);
s_t_b = sin(theta_b*rad);

% XL0=ATAN(H11/G11)
% SQQ=SQRT(G11**2+H11**2)
% SQ=SQRT(G11**2+H11**2+G10**2)
% ST0=SQQ/SQ
% CT0=G10/SQ
% SL0=-H11/SQQ
% CL0=-G11/SQQ
% CTCL= CT0*CL0
% CTSL=CT0*SL0
% STCL=ST0*CL0
% STSL=ST0*SL0

A = [[+c_t_b*c_p_b +c_t_b*s_p_b  -s_t_b]
    [    -s_p_b     +c_p_b     0]
    [+s_t_b*c_p_b +s_t_b*s_p_b  +c_t_b]];
if i_trans == -1; A = A'; end; % gm -> gg

c_t = cos(theta_gg*rad);
s_t = sin(theta_gg*rad);
c_p = cos(phi_gg*rad);
s_p = sin(phi_gg*rad);

z = c_t;
x = s_t .* c_p;
y = s_t .* s_p;

x_gm = A(1,1)*x + A(1,2)*y  +A(1,3)*z;
y_gm = A(2,1)*x + A(2,2)*y  +A(2,3)*z;
z_gm = A(3,1)*x + A(3,2)*y  +A(3,3)*z;

theta_gm = 90 - atan2(z_gm, sqrt(x_gm.^2 + y_gm.^2))/rad;
phi_gm = mod(atan2(y_gm, x_gm)/rad, 360);

if nargin == 5
    B_theta = varargin{2};
    B_phi   = varargin{3};
    BE = B_theta.*c_t;
    Bx = BE.*c_p - B_phi.*s_p;
    By = BE.*s_p + B_phi.*c_p;
    Bz = -B_theta.*s_t;
    
    Bx_gm = A(1,1)*Bx + A(1,2)*By  +A(1,3)*Bz;
    By_gm = A(2,1)*Bx + A(2,2)*By  +A(2,3)*Bz;
    Bz_gm = A(3,1)*Bx + A(3,2)*By  +A(3,3)*Bz;
    
    c_t = cos(theta_gm*rad);
    s_t = sin(theta_gm*rad);
    c_p = cos(phi_gm*rad);
    s_p = sin(phi_gm*rad);
    BE           = Bx_gm.*c_p + By_gm.*s_p;
%     B_r = BE.*s_t + Bz_gm.*c_t;
    varargout{1} = BE.*c_t    - Bz_gm.*s_t; % B_theta
    varargout{2} = By_gm.*c_p - Bx_gm.*s_p; % B_phi    

%     plot(B_r)
    
end
